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We present a new reciprocal space analytical method to cutoff the long range interactions in 
supercell calculations for systems that are infinite and periodic in 1 or 2 dimensions, extending 
previous works for finite systems. The proposed cutoffs are functions in Fourier space, that are 
used as a multiplicative factor to screen the bare Coulomb interaction. The functions are analytic 
everywhere but in a sub-domain of the Fourier space that depends on the periodic dimensionality. 
We show that the divergences that lead to the non-analytical behaviour can be exactly cancelled 
when both the ionic and the Hartree potential are properly screened. This technique is exact, fast, 
and very easy to implement in already existing supercell codes. To illustrate the performance of the 
new scheme, we apply it to the case of the Coulomb interaction in systems with reduced periodicity 
(as one- dimensional chains and layers) . For those test cases we address the impact of the cutoff in 
different relevant quantities for ground and excited state properties, namely: the convergence of the 
ground state properties, the static polarisability of the system, the quasiparticle corrections in the 
GW scheme and in the binding energy of the excitonic states in the Bethe-Salpeter equation. The 
results are very promising. 



PACS numbers: 02.70.-c,31.15.Ew,71.15.-m,71.15.Qe 

I. INTRODUCTION 

Plane waves expansions with periodic boundary con- 
ditions have been proven to be a very effective way 
to exploit the translational symmetry of infinite crystal 
solids, in order to calculate the properties of the bulk, by 
performing the simulations in one of its primitive cells 
only 1 . The use of plane waves is motivated by several 
facts. First, the translational symmetry of the poten- 
tials involved in the calculations is naturally and easily 
accounted for in reciprocal space, through the Fourier 
expansion. Second, very efficient and fast algorithms ex- 
ist (like FFTW^) that allow us to calculate the Fourier 
transforms very efficiently. Third the expansion in plane 
waves is exact, since they form a complete set, and it is 
only limited in practice by one parameter, namely, the 
maximum value of the momentum, that determines the 
size of the chosen set. Fourth, in many cases, the use 
of Born-von Karman periodic boundary conditions sup- 
plies a conceptually easy (though artificial) way to get rid 
of the dependence of the properties of a specific sample 
on its surface and shape, allowing us to concentrate on 
the bulk properties of the system in the thermodynamic 
limits 

However, mainly in the last decade, increasing inter- 
est has been developed in systems at the nano-scale, like 
tubes, wires, quantum-dots, biomolecules, etc., whose 
physical dimensionality is, for all practical purposes, less 
than three4> These systems are still 3-dimensional (3D), 
but their quantum properties are those of a confined sys- 
tem in one or more directions, and those of a periodic ex- 
tended system in the remaining directions. Other classes 



of systems with the same kind of reduced periodicity are 
the classes of the polymers, and of the solids with defects. 

Throughout this paper we call nD-periodic a 3D ob- 
ject, that can be considered infinite and periodic in n di- 
mensions, being finite in the remaining 3 — n dimensions. 
In order to simulate this kind of systems, a commonly 
adopted approach is the supercell approximation^ 

In the supercell approximation the physical system is 
treated as a fully 3D-periodic one, but a new unit cell 
(the supercell) is built in such a way that some extra 
empty space separates the periodic replica along the di- 
rection^) in which the system is to be considered as fi- 
nite. This method makes possible to retain all the advan- 
tages of plane waves expansions and periodic boundary 
conditions. Yet the use of a supercell to simulate objects 
that are not infinite and periodic in all the directions, 
leads to artifacts, even if a very large portion of vacuum 
is interposed between the replica of the system in the 
non-periodic dimensions. 

In fact, the straightforward application of the supercell 
method generates in any case fake images of the original 
system, that can mutually interact in several ways, af- 
fecting the results of the simulation. It is well known 
that the response function of an overall neutral solid of 
molecules is not equal, in general, to the response of the 
isolated molecule, and converges very slowly to it, when 
the amount of vacuum in the supercell is progressively 
increased AS 

For instance, the presence of higher order multipoles 
can make undesired images interact via the long range 
part of the Coulomb potential. In the dynamic regime, 
multipoles are always generated by the oscillations of the 
charge density even in systems whose unit cell does not 
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carry any multipole in its ground state. This is the case, 
for example, when we investigate the response of a system 
in presence of an external oscillating electric field. 

Things go worse when the unit cell carries a net charge, 
since the total charge of the infinite system represented 
by the supercell is actually infinite, while the charges at 
the surfaces of a finite, though very large system always 
generate a finite polarisation field. This situation is usu- 
ally normalised in the calculation by the introduction of 
a suitable compensating positive background charge. 

Another common situation in which the electrostat- 
ics is known to modify the ground state properties of 
the system occurs when a layered system is studied, and 
an infinite array of planes is considered instead of a sin- 
gle slab, being in fact equivalent to an effective chain of 
capacitors £ 

These issues become particularly evident in all the ap- 
proaches that imply the calculation of non-local opera- 
tors or response functions, because, in these cases, two 
supercells may effectively interact even if their charge 
densities do not overlap at all. This is the case, for exam- 
ple, of the many-body perturbation theory calculations 
(MBPT), and, in particular, of the self-energy calcula- 
tions at the GW level&i 

However we are usually still interested in the disper- 
sion relations of the elementary excitations of the system 
along its periodic directions, and those are ideally dealt 
with using a plane waves approach. Therefore, the ideal 
path to keep the advantages of the supercell formulation 
in plane waves, and to gain a description of systems with 
reduced periodicity free of spurious effects is to develop 
a technique to cut the Coulomb interaction off out of a 
desired region. This problem is not new and has been 
addressed now for a very long time and in different fields 
(condensed matter, classical fields, astrophysics, biology, 
etc). Several different approaches have been proposed in 
the past to solve it, however a complete review of them 
is beyond the scope of this paper. The aim of the present 
work is to focus on the widely used supercell schemes to 
show how the image interaction influences both the elec- 
tronic ground state properties and the dynamical screen- 
ing in the excited state of 0D-, 1D-, 2D-periodic systems, 
and to propose an exact method to avoid the undesired 
interaction of the replicas in the non-periodic directions. 

The paper is organised as follows: in Sec. [DJthe ba- 
sics of the plane wave method for solids are reviewed, in 
Sec. lIIII the new method is outlined, in Sec. lIVI the treat- 
ment of the singularities is explained, in Sec. [V] some 
applications of the proposed technique are discussed. 



II. THE 3D-PERIODIC CASE 

The main problem of electrostatics we are facing here 
can be reduced to the problem of finding solutions to the 
Poisson equation for a given charge distribution n(r), and 
given boundary conditions 

V 2 V(r) = -47rn(r). (1) 



In a finite system the potential is usually required to 
be zero at infinity. In a periodic system this condition 
is meaningless, since the system itself extends to infinity. 
Nevertheless the general solution of Eq. in both cases 
is known in the form of the convolution 



V(t) 



n(r') 



-dV 



(2) 



that it is referred from now on as the Hartree potential. 

It might seem that the most immediate way to build 
the solution potential for a given charge distribution is to 
compute the integral in real space, but problems immedi- 
ately arise for infinite systems. In fact, the density can be 
reduced to an infinite sum over delta charge distributions 
qS(r - r'), 



V(t) 



(3) 



and the integral in Eq. J2J becomes an infinite sum as 
well, but this sum is in general only conditionally, and 
not absolutely convergent^ The sum of Eq. © a poten- 
tial that is determined up to a constant for a neutral cell 
with zero dipole moment, while the corresponding sum 
for the electric field is absolutely convergent. A neutral 
cell with a non null dipole moment, on the opposite, gives 
a divergent potential, and an electric field that is deter- 
mined up to an unknown constant electric field (the sum 
for the electric field is conditionally convergent in this 
case). 

Even if, in principle, the surface terms have to be al- 
ways taken into account, in practice they are only rel- 
evant when we want to calculate energy differences be- 
tween states with different total charge. These terms can 
be neglected in the case of a neutral cell whose lowest 
nonzero multipole is quadrupoleai. As in the present 
work we are interested in macroscopic properties of the 
periodic system, those surface effects are never considered 
in the discussion that follows. However this sample-shape 
effects play an important role for the analysis of differ- 
ent spectroscopies as, for example, infrared and nuclear 
magnetic resonance. 

A major source of computational problems is the fact 
that the sum in Eq. J3J is very slowly converging when it 
is summed in real space, and this fact has historically mo- 
tivated the need for reciprocal space methods to calculate 
it. It was Ewald who first discovered that, by means of 
an integral transform, the sum can be split in two terms, 
and that if one is summed in real, while the other in 
reciprocal space, both of them are rapidly converging^ 
The point of splitting is determined by an arbitrary pa- 
rameter. 

Let us now focus on methods of calculating the sum in 
Eq. 10 purely based on the reciprocal space. 

If we consider a periodic distribution of charges with 
density n(r) such that n(r) = n(r + L n ), with L n = 
{n x li x ,n y Liy,n z 'L z }, and {n x ,n y ,n z } € Z, it turns out 
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that the reciprocal space expression for a potential like 



V{v) 



n(r')«(|r - r'|)d 



3,.' 



in a 3D-periodic system, can be written as 
V(G n ) = n(G n )u(G n ), 



(4) 



(5) 



where we have used the convolution theorem to trans- 
form the real space convolution of the density and the 
Coulomb potential into the product of their reciprocal 
space counterparts. Here G n = {n x G x ,n y G y ,n z G z } 
are the multiples of the primitive reciprocal space vec- 
tors j^-, and v(G n ) is the Fourier transform of the long 
range interaction v(r), evaluated at the point G n . For 
the Coulomb potential it is 

«(G») = %■ (6) 

Fourier transforming expression (|5j) back into real 
space we have, for a unit cell of volume f2, 



▼ r/ x 4vr ^ n(G„) 

V ^ = ~nl^ ~7p~ ex P0 G " • r )- 



At the singular point n x = n y — n z — the potential 
V is undefined, but, since the value at G = corresponds 
to the average value of V in real space, it can be chosen to 
be any number, corresponding to the arbitrariness in the 
choice of the static gauge (a constant) for the potential. 
Observe that the same expression can be adopted in the 
case of a charged unit cell, but this time, the arbitrary 
choice of v(G) in G = corresponds to the use of a 
uniform background neutralising charge. 



III. SYSTEMS WITH REDUCED PERIODICITY 

It has been shownii that the slab capacitance effect 
mentioned in the introduction actually is a problem that 
cannot be solved by just adding more vacuum to the su- 
percell. This has initially led to the development of cor- 
rections to Ewald's original methodic, and then to rig- 
orous extensions in 2D and lDi^i^ The basic idea is to 
restrict the sum in reciprocal space to the reciprocal vec- 
tors that actually correspond to the periodic directions 
of the system. These approach are in general of order 
0(A^ 2 )iLi^, but they have been recently refined to order 
0(N lniV)ii2i22i Another class of techniques, developed 
so far for finite systems, is based on the expansion of 
the interaction into a series of multipoles (fast multipolc 
method) pSiS&iSSi With this technique it is possible to eval- 
uate effective boundary conditions for the Poisson's equa- 
tions at the cell's boundary, so that the use of a supercell 
is not required at all, making it computationally very 
efficient for finit o 24 ' 25 and extended systems 2 ^. Other 



known methods, tipically used in molecular dynamics 
simulations, are the multipole-correction method2i, and 
the particle-mesh method^, whose review is beyond the 
scope of the present work, and we refer the reader to the 
original works for details. 

Differently from what happens for the Ewald sum, 
the method that we propose to evaluate the sum in 
Eq. @ entirely relies on the Fourier space and amounts 
to screening the unit cell from the undesired effect of 
(some of) its periodic images. The basic expression is 
Eq. (jSJl, whose accuracy is only limited by the max- 
imum value Gn of the reciprocal space vectors in the 
sum. Since there is no splitting between real and recip- 
rocal space, no convergence parameters are required. 

Our goal is to transform the 3D-periodic Fourier rep- 
resentation of the Hartree potential of Eq. JSJ into the 
modified one 



V(G n ) = n(G n )v(G n ) 



(8) 



such that all the interactions among the undesired peri- 
odic replica of the system disappear. The present method 
is a generalisation of the method proposed by Jarvis et 
alm& for the case of a finite system. 

In order to build this representation, we want to: 1) 
define a screening region V around each charge in the 
system, out of which there is no Coulomb interaction; 2) 
calculate the Fourier transform of the desired effective 
interaction v(r) that equals the Coulomb potential in D, 
and is outside V 



V(r) 



i if r € T> 
if r <£. V 



(9) 



Finally we must 3) modify the density n(r) in such a way 
that the effective density is still 3D-periodic, so that the 
convolution theorem can be still applied, but densities 
belonging to undesired images are not close enough to 
interact through v(r). 

The choice of the region T> for step 1) is suggested by 
symmetry considerations, and it is a sphere (or radius R) 
for finite systems, an infinite cylinder (of radius R) for 
ID-periodic systems, and an infinite slab (of thickness 
2R) for 2D-periodic systems. 

Step 2) means that we have to calculate the modified 
Fourier integral 



V(G) 



v(r)e~ 



iGr j3 



d^r = 



v(r)e 



iGr d 3 i 



spa ce 



( 10 ) 

Still we have to avoid that two neighbouring images in- 
teract by taking them far away enough from each other. 
Then step 3) means that we have to build a suitable su- 
percell, and re-define the density in it. 

Let us examine first step 2), i.e. the cutoff Coulomb 
interaction in reciprocal space. We know the expression 
of the potential when it is cutoff in a sphered It is 

4-7T 

\G) = -^[l-co S (GR)]. (11) 



v 0D < 
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The limit R — > oo converges to the bare Coulomb term in does not exist, since for G z = 0, the cutoff has a finite 



the sense of a distribution, while, since limc^o v 0D (G) = 
2nR 2 , there is no particular difficulty in the ori- 
gin. This scheme has been successfully used in many 
applicationa 10 i 23 i 25 i 27 i 28 . 

The ID-periodic case applies to systems with infinite 
extent in the x direction, and finite in the y and z di- 
rections. The effective Coulomb interaction is then de- 
fined in real space to be out of a cylinder of radius 
R having its axis parallel to the x direction. By per- 
forming the Fourier transformation we get the following 
expression for the cutoff coulomb potential in cylindrical 
coordinates >2i 



v 1D (G x ,G ± ) 



An 
G 2 



I + G±RJ 1 (G±R)K (\G X \R) 



\G x \RJo(GxR)Ki{\G x \R) , (12) 



where J and K arc the ordinary and modified cylindrical 
Bessel functions, and G± = y'G 2 + G 2 . 

It is easy to realise that, since the K functions damp 
the oscillations of the J functions very quickly, for all 
practical purposes this cutoff function only acts on the 
very first smaller values of G, while the unscreened ^ 
behaviour is almost unchanged for the larger values. 

Unfortunately, while the J n (Q functions have a con- 
stant value for £ = 0, and the whole cutoff is well defined 
for Gj_ = 0, the Kq(£) function diverges logarithmically 
for £ — > 0. Since, on the other hand, Ki(^) w for 
small £, 



v 1D 



(G x , Gj 



- \og(G x R) for G± >0,G X ^ 0+ 

(13) 

This means that the limit lim G ^ + V (G) does not exist 
for this cutoff function, and the whole G x = plane 
is ill-defined. We will come back to the treatment of the 
singularities in the next section. We notice that this loga- 
rithmic divergence is the common dependence one would 
get for the electrostatic potential of a uniformly charged 
ID wireSS. It is expected that bringing charge neutrality 
in place would cancel this divergence (see below). 

The 2D-periodic case, with finite extent in the z di- 
rection, is calculated in a similar manner. The effective 
Coulomb interaction is defined in real space to be out 
of a slab of thickness 2R symmetric with respect to the 
xy plane. In Cartesian coordinates we get 



v 2U (G h G z ) 
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G 2 
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-G„R 



|Gz 

Gii 



sm(\G z \R) 



\G Z \R) 



(14) 



where Gii 



G 2 



G 2 . 



In the limit R — > oo the unscreened potential ^ is 
recovered. Similarly to the case of ID, the limit G — > 



value, while it diverges in the limit G 







1 



v'' u (G h G z )~ for G,| > 0, G z -> 0+ 



(15) 



So far we haven't committed to a precise value of the 
cutoff length R. This value has to be chosen, for each 
dimensionality, in such a way that it avoids the interac- 
tion of any two neighbour images of the unit cell in the 
non-periodic dimension. 

In order to fix the values of R we must choose the size 
of the supercell. This leads us to the step 3) of our proce- 
dure. We recall that even once the long range interaction 
is cutoff out of some region around each component of the 
system, this is not sufficient yet to avoid the interaction 
among undesired images. The charge density has to be 
modified, or, equivalently, the supercell has to be built 
in such a way that two neighbouring densities along ev- 
ery non periodic direction do not interact via the cutoff 
interaction. 

It is easy to see how this could happen in the simple 
case of a 2D square cell of length L: if both r and r' 
belongs to the cell, then r, r 1 < L, and |r — r'| < \p2L (see 
the schematic drawing in Fig. ^| . If a supercell is built 
that is smaller than (1 + V%)L, there could be residual 
interaction, and the cutoff would no longer lead to the 
exact removal of the undesired interactions. 

Let us call Aq the unit cell of the system we are work- 
ing on, and A — {Ai, i — — oo, • • • , oo} the set of all the 
cells in the system. If the system is nD-periodic this set 
only includes the periodic images of Aq in the n periodic 
directions. Let us call B the set of all the non-physical 
images of the system, i.e. those in the non-periodic di- 
rections. Then A U B = R 3 . Obviously, if the system is 
3D-periodic A = K 3 , and B = 0. 

In general we want to allow the interaction of the elec- 
trons in Aq with the electrons in all the cells Ai € A, 
but not with those Bi £ B. To obtain this we define the 
supercell Go 3 Aq such that, Vi 



if r S A Q , and r' <E A t then |r - r'| e G 
if r e A Ql and r' e B % then |r - r'| ^ G 



(16) 



(see Fig. \l\ior a simplified 2D sketch). The new density 
n(r) is such that 



if r e Aq then n(r) = n(r), 

if r e Go, and r ^ Aq then n(r) = 0. 



(17) 



The size Lq of the super-cell in the non-periodic direc- 
tions depends on the periodic dimensionality of the sys- 
tem. In order to completely avoid any interaction, even 
in the case the density of the system is not zero at the 
cell border, it has to be 




for finite 

for ID-periodic 

for 2D-periodic 



(18) 
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(1+7(2))l' 



Figure 1: Schematic description for the supercell construc- 
tion in a 2D system. The upper sketch corresponds to the 
2D-periodic case (i.e. a 2D crystal). The middle sketch cor- 
responds to a OD-periodic system, and the bottom one to a 
ID-periodic. In the OD-periodic case the systems in different 
cells do not interact, while in the ID-periodic the chains do 
not interact, but within each chain the interaction of all of its 
elements is permitted.. 



the constants consistently for the two potentials. Since 
we know in advance that the sum must be finite, we can 
include, so to speak, all the infinities into these constants, 
provided that we find a method to separate out the long 
range part of both potentials on the same footing. 

In what follows we show how charge neutrality can be 
exploited to obtain the exact cancellations when oper- 
ating with the cutoff expression of Sec. II I II in Fourier 
space. 

The total potential of the system is built in the follow- 
ing way: we separate out first short and long range contri- 
butions to the ionic potential by adding and subtracting 
a Gaussian charge density n + (r) = Zexp(— a 2 r 2 ). The 
potential generated by this density is V + (r) — Z crf ^ . 
The ionic potential is then written as 



V(r) = AV(r) 



,erf (ar) 



(19) 



where a is chosen so that AV(r) is localised within a 
sphere of radius r a , smaller than the cell size. The ex- 
pression of the ionic potential in reciprocal space is 



V(G) = 4tt 



rsin(Gr) . . . j ex P 
v -AV(r)dr 



G 
To 1 



G 



G 2 



( 2 °) 

which, for G = gives a finite contribution from the first 
term, and a divergent contribution from the second term 



+ oo 

V(G = 0) = 4tt J r 2 AV(r)dr 



(21) 



Actually, since the required super-cell is quite large, a 
compromise between speed and accuracy can be achieved 
in the computation, using parallelepiped super-cell with 
Lc = 2L for all the cases. This approximation rests 
on the fact that the charge density is usually contained 
in a region smaller than the cell in the non-periodic di- 
rections, so that the spurious interactions are, in fact, 
avoided, even with a smaller cell. Therefore, on the basis 
of this approximation, we can choose the value of the cut- 
off length R always as half the smallest primitive vector 
in the non-periodic dimension. 



IV. CANCELLATION OF THE SINGULARITIES 

The main point in the procedure of eliminating the di- 
vergences in all the cases of interest is to observe that our 
final goal is usually not to obtain the expression of the 
Hartree potential alone, because all the physical quan- 
tities depend on the total potential, i.e. on the sum of 
the electronic and the ionic potential. When this sum is 
considered we can exploit the fact that each potential is 
defined up to an arbitrary additive constant, and choose 



The first is the contribution of the localised charge, and 
is easily computed, since the integrand is zero for r > r a . 
The second term is cancelled by the corresponding G = 
term in the electronic Hartree potential, due to the charge 
neutrality of the system. This trivially solves the problem 
of the divergences in 3D-periodic systems. 

Now let us consider a ID-periodic system. The Hartree 
term alone in real space is given by 



V(x,y,z)=J2 // n(G x ,y',z') 



x v(G x ,y -y',z- z')e iG * x ' dy'dz' . (22) 

Invoking the charge neutrality along the chain axis, we 
have that the difference between electron and ionic den- 
sities satisfies 



[nion{G x = 0, y, z) - n c \(G x = 0, y, z)) dydz = 0. 

(23) 

Unfortunately, the cutoff function in Eq. (|12fl is diver- 
gent for G x = 0. So the effective potential results in an 
undetermined • oo form. However, we can work out an 
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analytical expression for it by denning first a finite cylin- 
drical cutoff, but then bringing the size of the cylinder 
to infinity. This way, as a first step, we get a new cutoff 
interaction in a finite cylinder of radius R, and length h, 
assuming that h is much larger than the cell size in the 
periodic direction. In this case the modified finite cutoff 
potential includes a term 



v 1D 



{G x ,r, h) cx log 



h + Vh 2 + r 2 



(24) 



which, in turn gives, for the particular plane G x = 0, 



v 1D (G x = 0, G±) « -4tt / rJ (G ± r) log(r)dr 

o 



AnRlog(2h)^^. (25) 



The effective potential is now split into two terms, but 
only the second one depends on h. The second step is 
achieved by going to the limit h — » +oo, to obtain the 
exact infinite cutoff. By calculating this limit, we no- 
tice that only the second term in the right hand side 
of Eq. I|25|) diverges. This term is the one that can be 
dropped due to charge neutrality (in fact it has the same 
form for the ionic and electronic charge densities). Thus, 
for the cancellation to be effective in a practical imple- 
mentation, we have to treat on the same way both the 
ionic and Hartrce Coulomb contributions. Of course the 
first term in the right hand side Eq. I|25|l has always to 
be taken into account, affecting both the long and the 
short range part of the cutoff potentials. 

Following this procedure, we are able to get a consid- 
erable computational advantage, compared, e.g., to the 
method originally proposed by Spataru et al.^ since our 
cutoff is just an analytical function of the reciprocal space 
coordinates, and the evaluation of an integral for every 
value of G X ,G± is not needed. The cutoff proposed in 
Ref . l30l is actually a particular case of our cutoff, obtained 
by using the finite cylinder for all the components of the 
G vectors: in this case the quadrature in Eq. (|25|l has to 
be evaluated for each G Xl G y , and G z , and a convergence 
study in h is mandatory (see discussion in Sec. lVBl and 
Fig. 0). 

In the ID-periodic case, the G = value is now well 
defined, and it turns out to be \imG ± ~>o v(G x , G±) 



(G x = 0, G x = 0) = -7ri? 2 (2 log(fl) - 1). (26) 
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ID-periodic v 1D (G x , G±) = 
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+ G±RJi(G±R)K (G x R) 
-G X RJ (G ± R)K 1 (G X R)] 
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- 4tt J™ r.J {G±r) log(r)dr 








- 7ri? 2 (21og(i?) - 1) 


G n 


G z 


2D-periodic v {G\\,G Z ) = 
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any 


4^ [l + e - G li fl ( p. sin(G ZJ R) - cos{G z R))] 





> 


^ [1 - cos(G 2 J?) - G z Rsm(G z R)} 








-2ttR 2 



Table I: Reference Table summarising the results of the cut- 
off work for charge-neutral systems: finite systems (0D), one- 
dimensional systems (ID) and two-dimensional systems (2D). 
The complete reciprocal space expression of the Hartree po- 
tential is provided. For the ID case, R stands for the radius 
of the cylindrical cutoff whereas in the 0D case is the radius 
of the spherical cutoff. In 2D stands for half the thickness of 
the slab cutoff (see text for details). 



ratio 4^ between the in-plane lattice vectors. 



v 2U (G«=0,G z ) 



4tt 



I - cos(G z R) - G z Rsm(G z R)} 



8/1 log 



{a + VI + a 2 )(l + VI + a 2 ) \ sin(G z R) 



G z 



The G = value is 



v 2D 



(G,| =0,G Z = 0) 



~2nR 2 
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(28) 



To summarise, the divergences can be cancelled also 
in ID-periodic and 2D-periodic systems provided that 
I) we apply the cutoff function to both the ionic and 
the electronic potentials, 2) we separate out the infinite 
contribution as shown above, and 3) we properly account 
for the short range contributions as stated in Table [I] 
The analytical results of the present work are condensed 
in Tab. H] all the possible values for the cutoff functions 
are listed there as a quick reference for the reader. 



The analogous result for the 2D-periodic cutoff is ob- 
tained by imposing finite cutoff sizes h x — ah y = h 
(much larger than the cell size), in the periodic direc- 
tions x and y, and dropping the /i-dependent part before 
passing to the limit h — ► +00. The constant a is the 



V. RESULTS 

The scheme illustrated above has been implemented 
both in the real space time-dependent DFT code 



OCTOPUS , and in the plane wave many-body- 
perturbation-theory (MBPT) code SELF 32 . The tests 
have been performed on the prototypical cases of in- 
finite chains of atoms along the x axis. The compar- 
isons are performed between the 3D-periodic calculation 
(physically corresponding to a crystal of chains), and the 
ID-periodic case (corresponding to the isolated chain) 
both in the usual supercell approach, and within our ex- 
act screening method. The discussion for the 2D cases 
follows the same path as for the ID case, while results 
for the finite systems have already been reported in the 
literatureiSi 2 ^ We addressed different properties to see 
the impact of the cutoff at each level of calculation, from 
the ground state to excited state and quasiparticle dy- 
namics. 




y axis (a.u.) 



A. Ground state calculations 

All the calculations have been done with the real-space 
implementation of DFT in the OCTOPUS^ code. We 
have used non-local norm-conserving pseudopotentials 34 
to describe the electron-ion interaction and the local- 
density approximation (LDA) 35 to describe exchange- 
correlation effects. The particular choice of exchange- 
correlation or ionic-pseudopotential does not matter here 
as we want to assert the impact of the Coulomb cutoff 
and this is independent of those quantities. Moreover, 
we have used a grid of 0.38 a.u. for Si and Na. 

In this case the footprint of the interaction of neigh- 
bouring chains in the y and z direction is the dispersion 
of the bands in the corresponding direction of the Bril- 
louin zone. However it is known that, if the supercell is 
large enough, the bands along the T — X direction are un- 
changed. This is in apparent contradiction with the fact 
that the radial ionic potential for a wire (that asymp- 
totically goes like ln(r) as a function of the distance r 
from the axis of the wire) is completely different from 
the crystal potential. 

The answer to this contradiction is clear if we perform 
a cutoff calculation. In fact the overall effect on the oc- 
cupied states turns out to be cancelled by the Hartree 
potential, i.e. by the electron screening of the ionic po- 
tential, but two different scenarios are visible as soon as 
the proper cutoff is used. 

In Fig. |21 (top) it is shown the ionic potential, the 
Hartree potential, and their sum for a Si atom in a par- 
allelepiped supercell with side lengths of 2.5, 11, and 11 
a.u. respectively in the x, y and z directions. No cut- 
off is used here. The ionic potential is roughly behaving 
like i in the area not too close to the nucleus (where the 
pseudopotential takes over). The total potential, on the 
other hand, falls off rapidly to an almost constant value 
at around 4 a.u. from the nuclear position, by effect of 
the electron screening. 

Fig. |21 (bottom) shows the results when the cutoff is 
applied (the radius of the cylinder is R = 5.5 a.u. such 
that there is zero interaction between cells). The ionic 
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Figure 2: Calculated total and ionic and Hartree potentials 
for a 3D-periodic (top) and ID-periodic (bottom) Si chain 



potential now behaves like it is expected for a potential 
of a chain, i.e. diverges logarithmically, and is clearly 
different from the latter case. Nevertheless the sum of 
the ionic and Hartree potential is basically the same as 
for the 3D-periodic system. 

In the static case the two band structure are then ex- 
pected, and are found to be the same, confirming that, 
as far as static calculations are performed, the super- 
cell approximation is good, provided that the supercell 
is large enough (see Fig. 0). In static calculations, then, 
the use of our cutoff only has the effect of allowing us to 
eventually use a smaller supercell, what provides clear 
computational savings. In the case of the Si-chain a 
full 3D calculation would need of a cell size of 38 a.u. 
whereas the cutoff calculation would give the same result 
with a cell size of 19 a.u. Of course, when more delocal- 
ized states are considered, like higher energy unoccupied 
states, larger differences are observed with respect to the 
supercell calculation. 

In Fig. 0]a Na chain with lattice constant 7.5 a.u. is 
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Figure 3: Si linear chain in a supercell size of 4.9x19x19 a.u. 
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Figure 4: Effect of the cutoff in a Na linear chain in a su- 
percell size of 7.5x19x19 a.u. The bands obtained with an 
ordinary supercell calculation with no cutoff (dashed line) are 
compared to the bands obtained applying the ID cylindrical 
cutoff (solid line). As it is explained in the text, only the 
unoccupied levels are affected by the cutoff. 



considered in a cell of 7.5x19x19 a.u., and the effect of the 
cutoff on the occupied and unoccupied stated is shown. 
As expected, the occupied states are not affected by the 
use of the cutoff, since the density of the system within 
the cutoff radius is unchanged, and the corresponding 
band is the same as it is found for an ordinary 3D super- 
cell calculation with the same cell size. However there is 
a clear effect on the bands corresponding to unoccupied 
states, and the effect is larger the higher is the energy of 
the states. In fact the high energy states, and the states 
in the continuum are more delocalized, and therefore the 
effect of the boundary conditions is more sensible. 



B. Static polarisability 

After the successful analysis of the ground state prop- 
erties with the cutoff scheme, we have applied the mod- 
ified Coulomb potential to calculate the static polaris- 
ability of an infinite chain in the Random Phase Ap- 
proximation (RPA). As a test case we have considered 
a chain made of hydrogen atoms, two atoms per cell 
at a distance of 2 a.u. The lattice parameter was 4.5 
a.u. For this system we have also calculated excited state 
properties in many-body perturbation theory, in particu- 
lar the quasiparticle gap in Hcdin's GW approximation 6 
and the optical absorption spectra in the Bethe-Salpeter 
framework 7 - 36 (sec subsections below). All these calcu- 
lations have been performed in the code SELFi^ The 
polarisability for the monomer, i.e. a finite system, in 
the RPA approximation including local field effects is de- 
fined as 

a = -^o^ Xo0 ^- (29) 

where Xccil) ^ s the interacting polarisation function 
that is solution of the Dyson like equation 

xcG'(q) = xcG'(q) + Xl^GG"(q) u (q + G ")xc"G'(q)- 

G" 

(30) 

and x° is the non interacting polarisation function ob- 
tained by the Adler- Wiser expressions^ 7 . v(q+G) are the 
Fourier components of the Coulomb interaction. Note 
that the expression of a in Eq. |2Hlis also valid for calcu- 
lations in finite systems, in the supercell approximation, 
and the dependence from the wave-vector q is due to the 
representation in reciprocal space. 

In the top panel of Fig. we compare the values 
of the calculated polarisability a for different supercell 
sizes, a is calculated both using the bare Coulomb 
v(q + G) = q ^| 2 and the modified cutoff potential of 
Ea.lfT2j) (the radius of the cutoff is always set to half 
the inter-chain distance). The lattice constant along the 
chain axis is kept fixed. Using the cutoff the static polar- 
isability already converges to the asymptotic value with 
an inter-chain distance of 25 a.u., while without the cut- 
off the convergence is much slower, and the exact value is 
approximated to the same accuracy for much larger cell 
sizes (beyond the calculations shown in the top of Fig. El • 

We must stress that the treatment of the divergences 
in this case is different with respect to the case of the 
Hartree and ionic potential cancellation for ground-state 
calculations (i.e. charge neutrality). In fact, while in 
the calculation for the Hartree and ionic potential the 
diverging terms are simply dropped by virtue of the neu- 
tralising positive background, here the /i-dependence in 
Eq. H25|) can be removed only for the head component by 
virtue of the vanishing limit Hm q — >oXoo(q) = ^> wm l e f° r 
the other G x — components we have to resort to the 
expression of the finite cylindrical cutoff as in Ea. H25|) . 
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A finite version of the ID cutoff has been recently ap- 
plied to nanotube calculations. 30 ! 38 This cutoff was ob- 
tained by numerically truncating the Coulomb interac- 
tion along the axis of the nanotube, in addition to the 
radial truncation. Therefore the effective interaction is 
limited to a finite cylinder, whose size can be up to a 
hundred times the unit cell size, depending on the den- 
sity of the k-point sampling along the axis£ The cutoff 
axial length has to be larger than the expected bound 
exciton length. 

In the bottom part of Fig|5]we compare the results ob- 
tained with our analytical cutoff (Eq. l|12|l i with its finite 
special case as proposed in Ref. |3(J We observe that the 
value of the static polarisability calculated with the finite 
cutoff oscillates around an asymptotic value, for increas- 
ing axial cutoff lengths. The asymptotic value exactly 
coincides with the value that is obtained with our cutoff. 
We stress that we also resort to the finite form of the cut- 
off only for the diverging of components of the potential, 
thus we note that there is a clear numerical advantage in 
using our expression, since the cutoff is analytical for all 
values except at G x = 0, and the corresponding quadra- 
ture has to be numerically evaluated for these points only. 
In the inset of the bottom part of Fig. [S] it is also shown 
the convergences of the polarisability obtained with our 
cutoff with respect to the k-points sampling. The sam- 
pling is unidimensional along the axial direction. Observe 
that the calculation using our cutoff is already converged 
for a sampling of 20 k-points. In the upper axis it is also 
indicated the corresponding maximum allowed value of 
the finite cutoff length in the axial direction that has 
been used to calculate the G x = components. 

Finite size effects turn out to be relevant also for many- 
body perturbation theory calculations. For the same test 
system (linear H2-chain), in the next two subsections, 
we consider the performance of our cutoff potential for 
the calculation of the quasiparticle energies in GW— ap- 
proximation and in the absorption spectra in the Bcthc- 
Salpeter frameworki 3 ^ 



1. Quasiparticles in the GW approximation 

In the GW approximation, the non-local energy- 
dependent electronic self-energy £ plays a role sim- 
ilar to that of the exchange-correlation potential of 
DFT. T, is approximated by the convolution of the 
one electron Green's function and the dynamically 
screened Coulomb interaction W. We first calculate the 
ground state electronic properties using the DFT code 
ABINITi^ 3 . These calculation are performed in LDA 35 , 
and pseudopotentials 3 - 4 approximation. An energy cut- 
off of 30 hartree has been used to get converged re- 
sults. The LDA eigenvalues and eigenfunctions are then 
used to construct the RPA screened Coulomb interac- 
tion W, and the GW self-energy. The inverse dielec- 
tric matrix e^ G , has been calculated using the plasmon- 
pole approximation 39 and the quasiparticle energies have 
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Figure 5: Top: Polarisability per unit cell of an Hi chain 
in RPA approximation as a function of the supercell volume. 
The solid line joins the values obtained with the cutoff po- 
tential, while the dashed lines joins the values obtained with 
the bare Coulomb potential. The cutoff radius is 8.0 a.u. The 
inter-chain distance is indicated in the top axis. Bottom: Po- 
larisability of the H2 chain calculated with the finite cutoff 
potential of Refl.SOi In abscissa different values of the cutoff 
length along the chain axis. The dashed straight line indicates 
the value obtained with the cutoff of Ea, l)I2^ .rn the inset we 
show the convergence of the polarisability with respect to the 
k-points sampling along the chain axis obtained with the cut- 
off of Ea. ljf 2|l . In the upper axis it is indicated the maximum 
allowed length h for each k-point sampling used in the calcu- 
lation of the G x = components by Eq.JJSJ- 



been calculated at the first order of perturbation the- 
ory in £ — V xc m& Dividing the self-energy in an ex- 
change T, x and a correlation T, c parts ((4> F>FT \T,\(j)f FT ) — 

{<t>f FT \^f FT ) + (<t>f FT |S c |0f FT )), we get the follow- 
ing representation for the self-energy in a plane-waves 
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basis set: 



ik|S a (n,r 2 )|n'k') = -^ / t!§3E^+ G ) 

"Id. ' G 



Pnn 1 (q,G)p*,„ 1 (q,G)/„ lkl (31) 



and 



(nk|E c (r 1) r 2) c i ;)|nV) = ^ / £ w (q +G ') 

/d ' 



nl(k-q) 



1 _ /nl(k-q) 



r LDA 



uj — uj' — e"ZrV \ — iS u — uj' — e H, \ + i5 

nl(k— q) nl(k— q) 



(32) 



where p„ ni (q + G) = (nk|e l(q+G) ri |niki) and the in- 
tegral in the frequency domain in Ea. i|32[l has been an- 
alytically solved considering the dielectric matrix in the 



plasmon pole mode: ( e G G ,(w) 



J G,G' 



G,G 



~,2 



UJ 



G.G 



'))■ 



In order to eliminate the spurious interaction between 
different supercells, leaving the bare Coulomb interaction 
unchanged along the chain direction, we just introduce 
the expression of Eq. l|12[l in the construction of and 
S c , and also in the calculation of e 



GG' 



As we did for 

the calculation of the static polarisability, the divergences 
appearing in the components (G x = 0) cannot be fully 
removed and for such components we resort to the finite 
version of the cutoff potential Ea. H25|) . 

In Fig. [S] the convergence of the quasiparticle gap at 
the X point is calculated for different supercell sizes in 
the GW approximation. A cutoff radius of 8.0 a.u. has 
been used. When the cutoff potential is used, 60 k points 
in the axis direction has been necessary to get converged 
results. In the inset of Fig. E|we show the behaviour of 
the quasiparticle gap in function of the cutoff radius. We 
observe that for R c > 6 a.u. a plateau is reached, and, for 
R c > 12 a.u., a small oscillation appears due to interac- 
tion between the tails of the charge density of the system 
with its image in the neighbour cell. Differently from the 
DFT-LDA, calculation for neutral systems, where the su- 
percell approximation turns out to be good, as we have 
discussed above, we can see that the convergence of the 
GW quasiparticle correction turns out to be extremely 
slow with respect to the size of the supercell and huge 
supercells are needed in order to get converged results. 
This is due to the fact that in the GW calculation the 
addition of an electron (or a hole) to the system induces 
charge oscillation in the periodic images too. It is impor- 
tant to note that the slow convergence is caused by the 
correlation part of the self-energy ('Ea. (|32|l % ). while the ex- 
change part is rapidly convergent with respect to the cell 
size. The use of the cutoff Coulomb potential really im- 
proves drastically the convergence as it is evident from 
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Figure 6: Convergence of the GW quasiparticle gap for the 
H2 chain as a function of the cell size, using the bare Coulomb 
potential (dashed line) and the cutoff potential (solid line). In 
the inset the behaviour of the GW quasiparticle gap as a 
function of the value of the cutoff radius for a supercell with 
inter-chain distance of 32 a.u. is shown. The plateau obtained 
around a radius of 8 a.u. (i.e. one fourth of the supercell 
size) corresponds to the situation in which the radial images 
of chains no longer mutually interact, and the calculation is 
converged. Increasing the radius above approximately 12 a.u. 
the interaction is back and produces oscillations in the value 
of the gap. 



Fig© Notice that still at 38 a.u. inter-chain distance 
the GW gap is underestimated by about 0.5 eV. A sim- 
ilar trend (but with smaller variations) has been found 
by Onida et almM,, for a finite system (Sodium Tetramer) 
using the cutoff potential of Ea. (fTT|) . Clearly there is a 
strong dimensionality dependence of the self-energy cor- 
rection. The non-monotonic behaviour versus dimension- 
ality of the self-energy correction has also been pointed 
out in Ref. where the gap-correction was shown to 
have a strong component of the surface polarisation. 



2. Exciton binding energy: Bethe-Salpeter equation 

Starting from the quasiparticle energies we have calcu- 
lated the optical absorption spectra including electron- 
hole interactions solving the Bethe-Salpeter equation 36 . 
The basis set to describe the exciton state is composed by 
product states of the occupied and unoccupied LDA sin- 
gle particle states and the coupled electron-hole excited 
states \S) = ^2 cvk A cvk al k a vk \0) , where |0) is the ground 
state of the system. A cv -^ is the probability amplitude of 
finding an excited electron in the state (ck) and a hole 
in (uk), and it satisfies the equation 
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(33) 



Es is the excitation energy of the state \S) and K the 
interaction kernel that includes an unscreened exchange 
repulsive term K Exch and a screened electron-hole inter- 
action K dlr (direct term). In plane wave basis such terms 
read 
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(35) 



The screened potential has been treated in static RPA 
approximation (dynamical effects in the screening have 
been neglected as it is usually done in present Bcthc- 
Salpeter calculations 36 ). The quasiparticle energies en- 
tering in the diagonal part of the Hamiltonian in Eq. I|33l) 
are obtained applying a scissor operator to the LDA ener- 
gies, because in the studied test case the main difference 
between the quasiparticle and LDA band structure con- 
sists of a rigid energy shift of energy bands. From the 
solution of the Bethe-Salpeter equation (Ea l33|) it is pos- 
sible to calculate the macroscopic dielectric function, in 
particular the imaginary part reads 



£2(w) 
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where the summation runs over all the vertical excita- 
tions from the ground state |0) to the excited state \S), 
Es is the corresponding excitation energy, v is the veloc- 
ity operator and A is the polarisation vector. As in the 
case of GW calculation, in order to isolate the chain, we 
substitute the cutoff potential of Eci. (|12[) both in the ex- 
change term Eq. I|34l) and in the direct term of the Bcthc- 
Salpeter equation Ea. (|35|l . as well as for the RPA dielec- 
tric matrix present in Ea. (l35ll . 

In the top part of Fig. [JJwe show the calculated spectra 
for different cell sizes together with the non-interacting 
spectrum, and the spectrum obtained using the cutoff 
Coulomb potential for an inter-chain distance larger than 
20 a.u. and cylindrical cutoff radius of 8 a.u. The scissor 
operator applied in this calculation is the same for all the 
volumes and correspond to the converged GW gap. 

As it is known, the electron-hole interaction modifies 
both the shape and the energy of the main absorption 
peak. This effect is related to the the slow evolution 
of the polarisability per H2 uni1j2£. Furthermore, the 
present results clearly illustrate that the spectrum cal- 
culated without the cutoff slowly converges towards the 
exact result. This is highlighted in the bottom panel of 
Fig. [7J where we show the dependency of the exciton 
binding energy on the supercell volume, the binding en- 
ergy being defined as the energy difference between the 
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Figure 7: Top: Photo absorption cross section for different 
supercell volumes. In the legend the inter-chain distances 
corresponding to each volume are indicated. The intensity 
have been normalised to the volume of the supercell. The 
non-interacting absorption spectra and the spectra obtained 
with the cutoff potential are also included. Bottom: exciton 
binding energy vs supercell volume calculated using the cut off 
potential (solid line) , and the bare Coulomb potential (dashed 
line) . 



excitonic peak and the optical gap. We observe that the 
effect of the inter-chain interaction consists in reducing 
the binding energy with respect to its value in the iso- 
lated system. This value is slowly approached as the 
inter-chain distance increases, while, once the cutoff is 
applied to the Coulomb potential, the limit is reached 
as soon as the densities of the system and its periodical 
images do not interact. If we consider the convergence 
of the quasiparticle gap and of the binding energy with 
respect to the cell volume we notice that, if a cutoff is not 
used, the position of the absorption peak is controlled by 
the convergence of the Bethe-Salpeter equation solution, 
which, in turn, depends on the (slower) convergence of 
the GW energies. It is clear from Fig. that the use of 
the cutoff allows us to considerably speed up this bottle- 
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Figure 8: Values of e^, 1 ^, q y ) in the q z = plane for the Hi 
chain with an inter-chain distance of 25 a.u., using the bare 
Coulomb (top) and the cutoff potential (bottom). The axis 
of the chain is along the x direction 



fact allows us to sample the Brillouin zone only in the 
direction of the chain axis. For both the GW and Bcthc- 
Salpeter calculations a three dimensional sampling of the 
Brillouin zone is needed to get converged results when no 
cutoff is used, while a simple one-dimensional sampling 
can be adopted when the interaction is cutoff. 

VI. CONCLUSIONS 

An infinite system is an artifact that allows us to ex- 
ploit the powerful symmetry properties of an ideal system 
to approximate the properties of a finite one that is too 
large to be simulated at once. In order to use the valuable 
supercell approximation for systems that are periodic in 
less than three dimensions some cutoff technique is re- 
quired. The technique presented here provides a recipe 
to build the supercell, and a truncated Coulomb potential 
in Fourier space in such a way that the interactions to all 
the undesired images of the system are cancelled. This 
technique is exact for the given supercell sizes, and can be 
in most cases well approximated using supercells whose 
length is the double of the system size in the non-periodic 
dimensions. The method is very easily implemented in all 
available codes that use the supercell scheme, and is inde- 
pendent on the adopted basis set. We have tested it both 
in a real space code for LDA band-structure calculation 
of an atomic chain, and a plane wave code, for the static 
polarisability in RPA approximation, GW quasiparticle 
correction, and photo-absorption spectra in the Bethe- 
Salpeter scheme, showing that the convergence with re- 
spect the vacuum needed to isolate the system from its 
images is greatly enhanced, and the sampling of the Bril- 
louin zone is heavily reduced, being only necessary along 
the periodic directions of the system. 



neck. 

The use of our cutoff also has an important effect on 
the Brillouin zone sampling. In Fig. [HJwe show the value 
of e^Q (QxyQy) in the Qz = plane for a supercell corre- 
sponding to an inter-chain distance of 25 a.u. When the 
cutoff potential is used (bottom panel) the screening is 
smaller, compared to the case of the bare potential (top 
panel). Looking at the direction perpendicular to the 
chain (the chain axis is along the x direction) we see that 
the dielectric matrix is approximately constant, and this 
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